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Abstract 

We present the Unified Form Language (UFL), which is a domain-specific language for repre- 
senting weak formulations of partial differential equations with a view to numerical approx- 
imation. Features of UFL include support for variational forms and functionals, automatic 
differentiation of forms and expressions, arbitrary function space hierarchies for multi-field 
problems, general differential operators and flexible tensor algebra. With these features, UFL 
has been used to effortlessly express finite element methods for complex systems of partial 
differential equations in near-mathematical notation, resulting in compact, intuitive and read- 
able programs. We present in this work the language and its construction. An implementation 
of UFL is freely available as an open-source software library. The library generates abstract 
syntax tree representations of variational problems, which are used by other software libraries 
to generate concrete low-level implementations. Some application examples are presented and 
libraries that support UFL are highlighted. 

1 Introduction 

We present a language for expressing variational forms of partial differential equations (PDEs) in 
near-mathematical notation. The language, known as the Unified Form Language (UFL), inherits 
the typical mathematical operations that are performed on variational forms, thereby permitting 
compact and expressive computer input of mathematical problems. The complexity of the input 
syntax is comparable to the complexity of the classical mathematical presentation of the problem. 
The language is expressive in the sense that it provides basic, abstract building blocks which 
can be used to construct representations of complicated problems; it offers a mostly dimension- 
independent interface for defining differential equations; and it can be used to define problems 
that involve an arbitrary number of coupled fields. The language is developed with finite element 
methods in mind, but most of the design is not restricted to a specific numerical method. 

UFL is a language for expressing variational statements of partial differential equations and 
does not provide a problem solving environment. Instead, it generates abstract representations of 
problems that can be used by form compilers to create concrete code implementations in general 
programming languages. There exist a number of form compilers that generate low-level code 
from UFL. These include the FEniCS Form Compiler (FFC) [HI [HI EH E2], the SyFi Form 
Compiler (SFC) [7] and the Manycore Form Compiler [3(31 EI]- From a common UFL input, 
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these compilers differ in the strategies used to create and optimize a low-level implementation, 
and in the target low-level language. The code generated by these form compilers can be used in 
a problem solving environment, linked at compile time or dynamically at runtime. 

An example of a problem solving environment that uses code generated from UFL input is 
DOLFIN [511271, which is developed as part of the FEniCS Project [25]. Users of DOLFIN may 
describe a finite element discretization of a partial differential equation in UFL, and call a form 
compiler such as FFC or SFC to generate low-level code. In the case of FFC and SFC, this low- 
level code conforms to the UFC specification [El [2]) which is a C++ interface for functionality 
related to evaluation of local stiffness matrices, finite element basis functions and local-to-global 
mappings of degrees of freedom. The UFC code may then be used by DOLFIN to assemble and 
solve the linear or nonlinear system corresponding to the finite element discretization described in 
UFL. 

UFL is implemented as a domain-specific embedded language (DSEL) in Python. The dis- 
tinction between a DSEL and a high-level software component lies in the level of expressiveness; 
UFL expressions can be composed and combined in arbitrary ways within the language design 
limits. As an embedded language, UFL does not have a formal grammar of its own, but is instead 
implemented as a Python module and relies on the parser and grammar of the host language, 
Python. The UFL Python module defines types (classes) and operators that together form an 
expressive language for representing weak formulations of partial differential equations. In addi- 
tion, UFL provides a collection of algorithms for operating on UFL expressions. By implementing 
UFL as a DSEL in Python, we sacrifice some control over the syntax, but believe that this is 
overwhelmingly outweighed by the advantages. First, parsing is inherited and users may rely on 
all features of the Python programming language when writing UFL code, for example to define 
new operators. Second, it also permits the seamless integration of UFL into Python-based prob- 
lem solving environments. The Python interface of the library DOLFIN is an example of this. 
In particular, the use of just-in-time (JIT) compilation facilitates the incorporation of UFL in 
a scripted environment without compromising the performance of a compiled language. This is 
discussed in detail in Logg and Wells [24] . 

There have been a number of efforts to create domain-specific languages for scientific com- 
puting applications. Examples include SPL |£5] for signal processing and the Tensor Contraction 
Engine [14] for quantum chemistry applications. In the context of partial differential equations, 
there have been a number of efforts to combine symbolic computing, code generation and numer- 
ical methods. In some cases the code generation is explicit, while in other cases, such as when 
employing templates, implicit. Early examples include FINGER [15] and the Symbolic Mechanics 
System |22) . Analysa |13j is an abstract finite element framework of limited scope built upon 
Scheme. Feel++ [3H EH] uses C++ templates to create an embedded language for solving partial 
differential equations using finite element methods. Another example of a domain-specific lan- 
guage embedded in C++ is Sundance |28j . Sundance relies heavily on automatic differentiation 
to provide a problem solving environment targeted at PDE-constrained optimization. UFL also 
provides automated differentiation of functionals and variational forms, but the approach differs 
in some respects from Sundance. This is discussed later in this work. UFL is distinguished from 
the aforementioned efforts by its combination of a high level of expressiveness, mathematically- 
driven abstractions, extensibility, breadth of supported mathematical operations and embedment 
in a modern, widely- used and freely available language (Python). Moreover, it is deliberately 
decoupled from a code generator and problem solving environment. This provides modularity and 
scope to pursue different code generation and/or solution strategies from a common description 
of a variational problem. This is highlighted by the existence of the different form compilers that 
support UFL, with each targeting a specific code generation strategy or architecture. Unlike some 
of the efforts listed above, UFL is freely available under a GNU public license (LGPLv3+). 

The syntax used in UFL has its roots in FFC which was first released in 2005. At the time, FFC 
filled the roles of both form language and form compiler for the FEniCS Project. Much of the UFL 
syntax is inherited from early versions of FFC, but has since been re-implemented, generalized 
and extended to provide a more consistent mathematical environment, to cover a richer class of 
nonlinear forms and to provide a range of abstract algorithms, including differentiation. FFC no 
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longer provides an input syntax, rather it generates code from a UFL representation. The UFL 
form language was first released in 2009 [5] and has since then been tested on a wide range of 
applications. A rich and varied selection of applications that use UFL are presented in Logg et al. 

m- 

The remainder of this work is structured as follows. Section[2]summarizes the main mathemat- 
ical concepts on which UFL is based. A detailed presentation of the UFL language is then given 
in Section [3] This is followed in Section [4] by a number of examples that demonstrate the use of 
UFL for a variety of partial differential equations. The subsequent sections focus on the technical 
aspects of the UFL design. In Sections [5] and [6j we describe the internal representation of UFL 
expressions and provide an overview of the algorithms provided by UFL, respectively. Particular 
emphasis is placed on differentiation. Section [7] provides a brief discussion of validation and code 
correctness. Some conclusions are then drawn in Section [HI 

The implementation of UFL is available at https : //launchp ad.net/uf 1[ The e xamples pre- 
sented in this work, including the UFL code used, are archived at http://www.dspace.cam.ac. 
|uk/handle/1810/24398Tl 

2 Mathematical concepts and scope 

To clarify the notation, conventions, scope and assumptions of UFL and this paper, we begin 
by defining some key concepts in mathematical terms. We assume familiarity with variational 
formulations of PDEs and finite clement methods. These variational formulations are assumed 
to be expressed as sums of integrals over geometric domains. Each integrand is an expression 
composed from a set of valid functions and geometric quantities, with various operators applied. 
Each such function is an element of a function space, typically, but not necessarily, a finite element 
space, while the set of permitted operators include differential operators and operators from tensor 
algebra. The central mathematical abstractions, including multi-linear variational forms, tensor 
algebra conventions and the finite element construction, are formally introduced in the subsections 
below. 

When enumerating n objects, we count from 1 to n, inclusive, in the mathematical notation, 
while we count from to n — 1, inclusive, in computer code. 

2.1 Variational forms 

UFL is centered around expressing finite element variational forms, and in particular real-valued 
multi-linear forms. A real- valued multi-linear form a is a map from the product of a given sequence 
{Vj}j =1 of function spaces: 

a : V p x • • • x V 2 x Vi -> K, (1) 

that is linear in each argument. The spaces Vj are labeled argument spaces. For the case p < 2, 
V\ is referred to as the test space and V 2 as the trial space. The arity of a form p is the number of 
argument spaces. Forms with arity p = 0, 1, or 2 are named Junctionals, linear forms and bilinear 
forms, respectively. Such forms can be assembled on a finite element mesh to produce a scalar, 
a vector and a matrix, respectively. Note that the argument functions (vj)j are enumerated 
backwards such that their numbering matches the corresponding axis in the assembled tensor. 

If the form a is parametrized over one or more coefficient functions, we express the form as the 
mapping from a product of a sequence {VFfc}? =1 of coefficient spaces and the argument spaces: 

a : Wi x W 2 x • • • x W n x V p x • • • x V 2 X V x -> R, 
a>-^ a(w 1 ,w 2 ,...,w n ;v p ,...,v 2 ,v 1 ). 

Note that a is assumed to be (possibly) non-linear in the coefficient functions Wk and linear in the 
argument functions Vj . For a detailed exposition on finite element variational forms and assembly, 
we refer to [21] and references therein. To make matters concrete, we here list examples of some 
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a(u, v) := 


/ gradu • gradvda;, 
Jn 


p = 2, 


n = 


a(e; u, v) := 


/ e 2 gradu • gradvda;, 


p = 2, 


n = 1 


a(f;v) := 


/ fvdx, 
Jn 


p=l, 


n = 1 


a(u,v;) := 


/ grad(u — v) 2 da; 

7f2 


p = 0, 


n = 2 



forms with different arity p and number of coefficients n: 

(3) 
(4) 
(5) 
(6) 

where is the geometric domain of interest. 
2.1.1 Geometric domains and integrals 

UFL supports multi-linear forms defined via integration over geometric domains in the following 
manner. Let O C R d be a domain with boundary dil and let T = {T} be a suitable tessellation 
such that ft = {J Ter T. We denote the induced tessellation of dQ by T = {F}, and let J 70 denote 
the set of internal facets of T. Each of the three sets, the cells in T, the exterior facets in F and 
the interior facets in J 70 , is assumed to be partitioned into one or more disjoint subsets: 

n c n f n f 

T={J%, T={jT k , -F°=U^ (7) 

k=l fe=l fc=l 

where n c , n/ and n° denote the number of subsets of cells, exterior facets and interior facets, 
respectively. Given these definitions, it is assumed that the multi-linear form can be expressed in 
the following canonical form: 

a(w 1 ,w 2 ,...,w n ;v p ,...,v 2 ,v 1 ) = / Ik( w i> w 2,- ■ ■ ,w n ;v p ,. . . ,v 2 ,v 1 )dx 

+ / l[{wi,w 2 , ■ ■ ■ ,w n ;v p7 . . . ,v 2} v 1 )ds ,^ 

k=lF£F k JF 



n " 

+ / lt°(wi,w 2 ,...,w n ;v p ,...,v 2 ,v 1 )ds, 



where da; and ds denote appropriate measures. The integrand If. is integrated over the fcth subset 
7fc of cells, the integrand l[ is integrated over the kth subset of exterior facets and the integrand 
if' is integrated over the fcth subset J 7 ® of interior facets. 

UFL in its current form does not model geometrical domains, but allows integrals to be de- 
fined over subdomains associated with an integer index k. It is then the task of the user, or 
the problem-solving environment, to associate the integral defined over the subdomain k with a 
concrete representation of the geometrical subdomain. 



2.1.2 Differentiation of forms 

Differentiation of variational forms is useful in a number of contexts, such as the formulation of 
minimization problems and computing Jacobians of nonlinear forms. In UFL, the derivative of a 
form is based on the Gateaux derivative as detailed below. 

Let / and v be coefficient and argument functions, respectively, with compatible domain and 
range. Considering a functional M = M(/), the Gateaux derivative of M with respect to / in the 
direction v is defined by 

M'(f-v) ee D f M(f)[v] = ±- [M(f + rv)] T=0 . (9) 
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Given a linear form L( f; v) (which could be the result of the above derivation) and another 
compatible argument function it, we can continue by computing the bilinear form -£/(/; it, v); that 
is, the derivative of L with respect to / in the direction u, defined by 



L'(f;u,v) = D f L(f;v)[u} = — [L(f ■ 



C=o- 



(10) 



In general, this process can be applied to forms of general arity p > to produce forms of 
arity p + 1. Note that if the form to be differentiated involves an integral, we assume that the 
integration domain does not depend on the differentiation variable. To express the differentiation 
of a general form, consider the following compact representation of the canonical form Q: 



u JD k 



dpk, 



(11) 



where {-Dfc} and {dpk} are the geometric domains and corresponding integration measures, and 
{Ik} are the integrand expressions. We can then write the derivative of the general form (111 with 
respect to, for instance, w% in the direction v p+ i as 



D Wl F({w i )? =1 ;{v j ) 1 j=P )[v P+ i]=Y, 



dT I 



h(wi +tv p+ i, (tUi)" =2 ; (vj)j =p ) _ dp k . (12) 



T = 



2.2 Tensors and tensor algebra 

A core feature of UFL is its tensor algebra support. We summarize here some elementary tensor 
algebra definitions, notation and operations that will be used throughout this paper. 

First, an index is either a fixed positiv^ integer, in which case it is labeled a fixed-index, or 
a symbolic index with no value assigned, in which case it is called a free-index. A multi-index 
is an ordered tuple of indices, each of which can be free or fixed. Moreover, a dimension is a 
positive fixed-index. A shape s is an ordered tuple of zero or more dimensions: s = (s l7 . . . , s r ); 
the corresponding rank r equals the length of the shape tuple. 

Any tensor can be represented either as a mathematical object with a (tensor) shape or in terms 
of its scalar components with reference to a given basis. More precisely, the following notation 
and bases are used. Scalars are considered rank zero tensors. We denote by {e z }f =1 the standard 
orthonormal Euclidean basis for ]R d of dimension d. A basis {E a } a for M s , where s = (si, . . . , s r ) 
is a shape, is naturally defined via outer products of the vector basis: 

E a = e ai ® • • • ® e Q " , (13) 

where the range of the multi-index a — (a%, . . . ,a r ) is such that 1 a, ^ s, for i — 1, . . . ,r. 
In general, whenever a multi-index a is used to index a tensor of shape s, it is assumed that 
1 on Si for i = 1, . . . , r. Then, the scalar component of index i of a vector v defined relative 
to the basis {e z }i is denoted Uj. More generally, for a tensor C of shape s, C a denotes its scalar 
component of multi-index a with respect to the basis {E a } a of M s . Moreover, whenever we 
write we i m ply ^2i—i v i: where d is the dimension of v. Correspondingly, for sums over 

multi-indices, ^2 a C a implies Y] a • • ■ Y] a C a with the deduced ranges. 

Whenever one or more free-indices appear twice in a monomial term, summation over the 
free-indices is implied. Tensors v, A and C of rank 1, 2 and r, respectively, can be expressed using 
the summation convention as: 

v = Vl e\ A — AijE i] , C = C a E a . (14) 

We will also consider general tensor- valued functions / : il —> M s , 51 C M. d , where the shape 
s in this context is termed the value shape. Indexing of tensor-valued functions follows the same 

1 Indices are positive in the mathematical base 1 notation used here and non-negative in the base notation 
used in the computer code. 
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notation and assumptions as for tensors. Furthermore, derivatives with respect to spatial coordi- 
nates may be compactly expressed in index notation with the comma convention in subscripts. For 
example, for coordinates (x±, . . . , Xd) G fi and a function / : fi — ^ M, a vector function v : SI — > R" 
or a tensor function C : Jl — ¥ ]R S , we write for indices i and j, and multi- indices a and /3, with the 
length of j3 denoted by p: 

dxi ' 4J ctej ' dxp 1 ■ ■ ■ dxp p 

2.3 Finite element functions and spaces 

A finite element space Vh is a linear space of piecewise polynomial fields defined relative to a 
tessellation Th = {T} of a domain Q. Such spaces are typically defined locally; that is, each field 
in the space is defined by its restriction to each cell of the tessellation. More precisely, for a finite 
element space Vh of tensor-valued functions of value shape s, we assume that 

V h = {v G H : v\ T G V t }, V t = {v : T -> R s : w Q G P(T) Va}, (16) 

where the space "H indicates the global regularity and V — V(T) is a specified (sub-)space of 
polynomials of degree q ^ defined over T. In other words, the global finite element spaces Vh is 
defined by patching together local finite element spaces Vt over the tessellation Th- Note that the 
polynomial spaces may vary over the tessellation; however, this dependency is usually omitted for 
the sake of notational brevity. 

The above definition may be extended to mixed finite element spaces. For given local finite 
element spaces {Vj}™ =1 of respective value shapes {s l }™ =1 , we define the mixed local finite element 
space W by: 

W = Vi x V 2 x ••• x V„ = {w = (vi,v 2 ,...,v n ) : u< G V u i = 1, 2, . . . ,n}. (17) 



The extension to global mixed finite element spaces follows as in (16 1. Note that all element 
factors in a mixed element are assumed to be defined over the same cell T. The generalization to 
nested hierarchies of mixed finite elements follows immediately, and as such are also admitted. 

Any w G W has the representation w : T —> M*, with a suitable shape tuple t, and the value 
components of w must be mapped to components of v\ G Vj. Let r l be the corresponding rank 
of the value shape s*, and denote by = n^=i s j ^ ne corresponding waZue size; i.e, the number 
of scalar components. In the general case, we choose a rank 1 shape t = (X^iP*)' an< ^ ma P the 
flattened components of each v 1 to components of the vector-valued w; that is, 

w = (v 1 1 ,...,v 1 pl ,...,v^...,v; n ). (18) 

This mapping permits arbitrary combinations of value shapes s l . In the case where Vi coincide for 
all i = 1, . . . , n, we refer to the resulting specialized mixed element as a vector element of dimension 
n, and choose t = (n, s\, . . . , s^). As a generalization of vector elements, we allow tensor elements^] 
of shape c with rank q, which gives the value shape t = (ci, . . . , c q , s\, . . . , sh). Tensor elements 
built from scalar subelements may have symmetries, represented by a symmetry mapping from 
component to component. The number of subelements then equals n = (fliLi c i) ~ m > where m 
is the number of value components that are mapped to another. 

For some finite element discretizations, it is helpful to represent the local approximation space 
as the enrichment of one element with another. More precisely, for local finite element spaces 
Vi, V2, ■ ■ • , V n defined over a common cell T and of common value shape s, we define the space 

W = Vi+V 2 + ---+V n = {vi + v 2 + --- + v n : Vi G V, i = l,...,n}. (19) 



Again, the extension to global enriched finite clement spaces follows as in (16 1. The MINI ele- 
ment [11) for the Stokes equations is an example of an enriched element. 



2 Tensor- valued elements, not to be confused with tensor product elements. 
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Finite clement class specification 



FiniteElement (family, cell, degree) 
VectorElement (family, cell, degree, (dim)) 
TensorElement (family, cell, degree, (shape), (symmetry)) 



MixedElement (elements) 
EnrichedElement (elements) 
RestrictedElement (element , domain) 



Table 1: Overview of 
element classes defined 
in UFL. In addition 
to the arguments given 
here, a specific quadra- 
ture scheme can be given 
for the primitive finite 
elements (those defined 
by family name). Argu- 
ments in parentheses are 
optional. 



3 Overview of the language 

UFL can be partitioned into sublanguages for finite elements, expressions, and forms. We will 
address each separately below. Overall, UFL has a declarative nature similar to functional pro- 
gramming languages. Side effects, sequences of statements, subroutines and explicit loops found 
in imperative programming languages are all absent from UFL. The only branching instructions 
are inline conditional expressions, which will be further detailed in Section |3.2.6 



3.1 Finite elements 

The UFL finite element sublanguage provides syntax for finite elements and operations over finite 
elements, including mixed and enriched finite elements, as established in Section [273} 

3.1.1 Finite element abstractions and classes 

UFL provides four main finite clement abstractions: primitive finite elements, mixed finite ele- 
ments, enriched finite elements and restricted finite elements. Each of these abstractions provides 
information on the value shape, the cell and the embedding polynomial degree of the element 



(see (16)), and each is further detailed below. We remark that UFL is primarily concerned with 
properties of local finite element spaces: the global continuity requirement and the specific imple- 
mentation of the element degrees of freedom or basis functions are not covered by UFL. For an 
overview of finite element abstractions with initialization arguments, see Table [T] Example usage 
will be presented in Section [4] 

In the literature, it is common to refer to finite elements by their family parametrized by cell 
type and order: for instance, the Nedelec face elements of the second kind over tetrahedra of second 
order |34j . The global continuity requirements are typically implicitly implied by the family: for 
instance, it is generally assumed that the aforementioned Nedelec face element functions indeed 
do have continuous normal components across faces. Moreover, finite elements may be known by 
different family names, for instance the aforementioned Nedelec face elements coincide with the 
Brezzi-Douglas-Marini elements on tetrahedra [IH] , which again coincide with the "PA 2 (T) family 
on tetrahedra [12] , 

UFL mimics the literature in the sense that primitive finite elements are defined in terms of a 
family, a cell and a polynomial degree via the FiniteElement class (see Table [T]). Additionally, a 
quadrature scheme label can be given as an optional argument. The family must be an identifying 
string, while the cell is a geometric description of the domain. The UFL documentation contains 
the comprehensive list of preregistered families and cells. Multiple names in the literature for the 
same finite element are handled via family aliases. UFL supports finite element exterior calculus 
notation for simplices in one, two or three dimensions via such aliases. By convention, elements 
of a finite element family are numbered in terms of their polynomial degree q such that their 
fields are indeed included in the complete polynomial space of degree q. This facilitates internal 
consistency, although it might conflict with some notation in the literature. For instance, the 
lowest order Raviart-Thomas elements have degree 1 in UFL. 
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Operation 



Equivalent syntax 



M 


= U * V 


M 


= MixedElement(U, V) 


M 


= U * V * w 


M 


= MixedElement((U, V), W) 


M 


= U + V 


M 


= EnrichedElement (U, V) 


M 


= V['facet'] 


M 


= RestrictedElement(V, 'facet') 



Table 2: An overview of UFL op- 
erators over elements: examples 
of operator usage matched with 
the equivalent verbose syntax. 



Syntax is provided for defining vector elements. The VectorElement class accepts a family, a 
cell, a degree and the dimension of the vector element. The dimension defaults to the geometric 
dimension d of the cell. The value shape of a vector element is then (d, s) where s is the value shape 
of the corresponding finite element of the same family. Moreover, further structure can be imposed 
for (higher-dimensional) vector elements with a rank two tensor structure. The Tens or Element 
class accepts a family, a cell and a degree, and in addition a shape and a symmetry argument. 
The shape argument defaults to the tuple (d, d) and the value shape of the tensor element is 
then (d, d, s). The symmetry argument may be boolean true to define the symmetry Ay = Aji if 
the value rank is two. It may also be a mapping between the component tuples that should be 
equal, such as {(0,0) : (0,1) , (1,0) : (1,1)} to define the symmetries An — Ai 2 , A 2 i = A 22 - 
The vector and tensor element classes can be viewed as optimized, special cases of mixed finite 
elements. 

In general, mixed finite elements in UFL are created from a tuple of subelements through 
the MixedElement class. Each subelement can be a finite, vector, or tensor element as described 
above, or in turn a general mixed element. The latter can lead to nested mixed finite elements of 
arbitrary, though finite, depth. All subelements must be defined over the same geometric cell and 
utilize the same quadrature scheme (if prescribed) . The degree of a mixed finite element is defined 
to be the maximal degree of the subelements. Note that mixed finite elements are recursively 
flattened. Their value shape is (s, ) where s is the total number of scalar components. 

Enriched elements can be defined via the EnrichedElement class, given a tuple of finite, vector, 
tensor, or mixed subelements. The subelements must be defined on the same cell and have the 
same value shape. These then define the cell and value shape of the enriched element. The degree 
is inferred as the maximal degree of the subelements. 

Finally, UFL also offers a restricted element abstraction via the RestrictedElement class, 
taking as arguments any of the element classes described above and a domain. The term restricted 
in this setting refers to the elimination of element functions that vanish on the given domain; 
Labeur and Wells [23] provide an example utilizing elements restricted to cell facets. The value 
shape, cell and degree of a restricted element are directly deduced from the defining element. 

3.1.2 Operators over finite elements 

For readability and to reflect mathematical notation, UFL provides some operators over the finite 
element classes defined in the previous section. These operators include the binary operators 
multiplication (*) and addition (+), and an indexing operator ([]). These operators and their 
long-hand equivalent are presented in Table [2j The multiplication operator acts on two elements 
to produce a mixed element with the two elements as subelements in the given order. Note that the 
multiplication operator (in Python) is binary, so multiplication of three or more elements produces 
a nested mixed element. Similarly, the addition operator acts on two elements to yield an enriched 
element with the two given elements as subelements. Finally, the indexing operator restricts an 
element to the domain given by the index argument, thus returning a restricted element. 

3.2 Expressions 

The language for declaring expressions consists of a set of terminal expression types and a set of 
operators acting on other expressions. Each expression is represented by an object of a subclass of 
Expr. Each operator acts on one or more expressions and produces a new expression. An operator 
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result is uniquely defined by its operator type and its operand expressions, and cannot have non- 
expression data associated with it. A terminal expression does not depend on any other expressions 
and typically has non-expression data associated with it, such as a finite element, geometry data or 
the values of literal constants. Terminal expression types are subclasses of Terminal and operator 
results are represented by subclasses of Operator, both of which are subclasses of Expr. Any UFL 
expression object is the root of a self-contained expression tree in which each tree node is an Expr 
object. The references from objects of Operator subtypes to the operand expressions represent 
directed edges in the tree and objects of Terminal subtypes terminate the tree. 

As an embedded language, UFL allows the use of Python variables to store subexpression refer- 
ences for reuse. However, UFL itself does not have the concept of mutable variables. In fact, a key 
property of all UFL expressions, including terminal types, is their immutable statej^] Immutable 
state is a prerequisite for the reuse of subexpression objects in expression trees by reference instead 
of by copying. This aspect is critical for an efficient symbolic software implementation. 

The dependency set of an expression is the set of non-literal terminal expressions that can be 
reached from the expression root. An expression with an empty dependency set can be evaluated 
symbolically, but in general the evaluation of a UFL expression can only be carried out when the 
values of its dependencies are known. Numerical evaluation of the symbolic expression without 
code generation is possible when such values are provided, but this is an expensive operation and 
not suitable for large scale numerical computations. 

Every expression is considered to be tensor-valued and its shape must always be defined. Fur- 
thermore, every expression has a set of free-indices. Note that the free-index set of any particular 
expression object is not associated with its shape; for instance, if A is a rank two tensor with shape 
(3, 3) (and no free indices), then Ay is a rank zero tensor expression; in other words, scalar- valued 
and with the associated free-indices i and j. Mathematically one could see A and Ay as being the 
same, but represented as objects in software they are distinct. While A represents a matrix-valued 
expression, Ay represents any scalar values of A. Because the tensor properties of all subexpres- 
sions are known, dimension errors and inconsistent use of free-indices can be detected early. The 
following sections describe terminal expressions, the index notation and various operators in more 
detail. 

3.2.1 Terminal expressions 

Terminal expressions in UFL include literal constants, geometric quantities and functions. In 
particular, UFL provides a domain-specific set of types within these three fairly generic groups. 
Table [3] provides an overview of the literal constants, geometric quantities and functions available. 

Literal constants include integer and real- valued constants, the identity matrix, the Levi-Civita 
permutation symbol and unit tensors. Geometric quantities include spatial coordinates and cell 
derived quantities, such as the facet normal, facet area, cell volume, cell surface area and cell 
circumradius. Some of these are only well defined when restricted to facets, so appropriate errors 
are emitted if used elsewhere. 

Functions are cell- wise or spatially varying expressions. These are central to the flexibility of 
UFL. However, in contrast to other UFL expressions, functions are merely symbols or placeholders. 
Their values must generally be determined outside of UFL. All functions are defined over function 
spaces, introduced in Section [3.1[ such that their tensor properties, including their shape, can be 
derived from the function space. Functions are further grouped into coefficient functions and argu- 
ment functions. Expressions must depend linearly on any argument functions; see Section [2. 1| No 
such limitations apply to dependencies on geometric quantities or coefficient functions. Functions 
are counted, or assigned a count, as they are constructed, and so the order of construction matters. 
In particular, different functions are assumed to have different counts. For convenience, UFL pro- 
vides constructors for argument functions called TestFunction and TrialFunction which apply 
a fixed ordering. 

3 In the PyDOLFIN library, UFL function types are subclassed to carry additional mutable state which does not 
affect their symbolic meaning. UFL algorithms carefully preserves this information. 
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Literal tensor constants 



Mathematical notation UFL notation 

/ I = Identity (2) 

e eps = PermutationSymbol (3) 

e x ,e y ex, ey = unit_vectors (2) 

e x £x, SxSy, e y e x , e y e y exx, exy, eyx, eyy = unit_matrices (2) 



Geometric quantities Functions 



Math. 


UFL notation 


Math. 


UFL notation 


X 


x = cell.x 


ceK 


c = Constant (cell) 


n 


n = cell.n 


ge V 


g = Coefficient (element) 


\T\ 


h = cell. volume 


w e V 


w = Argument (element) 


r(T) 


r = cell . circumradius 


ueV 


u = TrialFunction(element) 


\F\ 


fa = cell . f acetarea 


v e V 


v = TestFunction(element) 


Efct\F\ 


ca = cell . cellsurf acearea 







Tabic 3: Tables of terminal expressions. The examples are given with reference to a predefined cell 
T C R d , with coordinates x € R d , facets {F} and facet normal n, denoted cell, and a predefined 
local finite element space V of some finite element denoted element. | • | denotes the volume, while 
r(T) denotes the circumradius of the cell T; that is, the radius of the circumscribed sphere of the 
cell. 



In the FEniCS pipeline, functions are evaluated as part of the form compilation or the assembly 
process. Argument functions are interpreted as a representation of each finite element basis func- 
tion. Moreover, the ordering of the argument functions defined by their count determines which 
global tensor axis each argument is associated with when assembling the global tensor (such as a 
sparse matrix) from a form. Form compilers may specialize the evaluation of the basis functions 
during compilation, or keep the choice of element space open until runtime. On the other hand, co- 
efficient functions are used to represent global constants, finite element fields or any function that 
can be evaluated at spatial coordinates during finite element assembly. The limitation to functions 
of spatial coordinates is necessary for the integration of forms to be a cell-wise operation. 

3.2.2 Index notation 

UFL mirrors conventional index notation by providing syntax for defining fixed and free-indices 
and for indexing objects. For convenience, free-index objects with the names i, j, k, 1, p, 
q, r, s are predefined. However, these can be redefined and new ones created. A single free- 
index object is created with i = Index (), while multiple indices are created with j, k, 1 = 
indices (3) . 

The main indexing functionality of UFL is summarized in Table [4j The indexing operator [] 
applied to an expression yields a component-wise representation. For instance, for a rank two 
tensor A (A) and free-indices i, j, A[i, j] yields the component- wise representation Aij. The 
mapping from a scalar-valued expression with components identified by free-indices to a tensor- 
valued expression is performed via as_tensor(A[i, j] , (i, j)). The as_vector, as_matrix, 
as_tensor functions can also be used to construct tensor-valued expressions from explicit tuples 
of scalar components. Note how the combination of indexing and as_tensor allows the reordering 
of the tensor axes in the expression Ayjy = B^uj. Finally, we remark that fixed and free-indices 
can be mixed during indexing and that standard slicing notation is available. 
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Mathematical notation 



UFL notation 



Bijki 
(a, b, c) 

A = Bm, (A, = Bi) 

A = BjiE ij , (Ai^Bji) 
iki 



A[i] 

B[i,j,k,l] 

as_vector ( (a, b, c)) 
as_vector (B [i] , i) 
as_matrix(B [j , i] , (i,j)) 



A = B kHj E^ ki , (Aijki = B kUj ) as_tensor(B[k,l,i, j] , (i,j,k,l)) 



Table 4: Table of indexing operators: i, j, k, 1 are free- indices, while a, b, c are other expressions. 



Mathematical notation UFL notation 



A + B 
A ■ B 
A : B 

AB = A ( 
Ax B 



A T 

sym A 
skew A 
dev A 
tr = At 



B 



A + B 
dot (A, B) 
inner (A, B) 
outer (A, B) 
cross(A, B) 



transpose (A) 
sym(A) 
skew (A) 
dev (A) 
tr(A) 



Mathematical notation 



det A 
cofac A 
A- 1 



v 
A 



Vi = An (no sum) 



A, 



An — 



Bij, if i = j, 
0, otherwise 
v h if i =j, 
0, otherwise 



UFL notation 



det (A) 
cofac (A) 
inv(A) 



diag_vector (A) 
diag(B) 

diag(v) 



Table 5: Table of tensor algebraic operators 



3.2.3 Arithmetic and tensor algebraic operators 

UFL defines arithmetic operators such as addition, multiplication, I 2 inner, outer and cross prod- 
ucts and tensor algebraic operators such as the transpose, the determinant and the inverse. An 
overview of common operators of this kind is presented in Table [5j These operators will be familiar 
to many readers and we therefore make only a few comments below. 

Addition and subtraction require that the left and right operands have the same shape and 
the same set of free-indices. The result inherits those properties. 

In the context of tensor algebra, the concept of a product is heavily overloaded. Therefore, 
the product operator * has no unique intuitive definition. Our choice in UFL is to define * as the 
product of two scalar- valued operands, one scalar and one tensor valued operand, and additionally 
as the matrix- vector and matrix-matrix product. The inner function defines the inner product 
between tensors of the same shape, while dot acts on two tensors by contracting the last axis of 
the first argument and the first axis of the second argument. If both arguments are vector- valued, 
the action of dot coincides with that of inner. 

When applying the product operator * to operands with free indices, summation over repeated 
free-indices is implied. Implicit summation is only allowed when at least one of the operands is 
scalar- valued, and the tensor algebraic operators assume that their operands have no repeated 
free- indices. 

3.2.4 Nonlinear scalar functions 

UFL provides a number of familiar nonlinear scalar real functions, listed in Tables [6] and [7| All 
of these elementary, trigonometric or special functions assume a scalar-valued expression with no 
free- indices as argument. Their mathematical meaning is well established and implementations 
are available in standard C++ or Boost [Tj. To apply any scalar function to the components of a 
tensor- valued expression, the element-wise operators listed in Table [7] can be used. 
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Mathematical notation UFL notation 



Mathematical notation UFL notation 



a/b 


a/b 


cos/ 


cos (f ) 


a b 


a**b, pow(a,b) 


sin/ 


sin(f) 


Vf 


sqrt(f ) 


tan/ 


tan(f ) 


exp/ 


exp(f ) 


arccos / 


acos (f ) 


In/ 


ln(f) 


arcsin / 


asin(f ) 


l/l 


abs(f ) 


arctan / 


atan(f ) 


sign/ 


sign(f ) 







Table 6: (Left) Table of elementary, nonlinear functions. (Right) Table of trigonometric functions. 



Mathematical notation UFL notation Mathematical notation UFL notation 



erf / 


erf (f ) 




C 


Cq 


— A a B a 


elenumilt (A , B) 


JAf) 


bessel_J(nu, 


f) 


C 


Cq 


= A a j B a 


elem_div(A, B) 


W) 


bessel_Y(mi, 


f) 


c 


C a 


- A B <* 

— A o> 


elem_pow(A, B) 


W) 


bessel_I (nu, 


f) 


c 


C a 


= f(A a ,...) 


elem_op(f , A, . . . ) 


KM) 


bessel_K(mi, 


f) 











Table 7: (Left) Table of special functions. (Right) Table of element-wise operators. 



3.2.5 Differential operators and explicit variables 

UFL supports a range of differential operators. Most of these mirror common ways of expressing 
spatial derivatives in partial differential equations. A summary is presented in Table [8] (left) . 

Basic partial derivatives, d/dxi, can be written as either A.dx(i) or Dx(A,i), where i is 
either a free-index or a fixed-index in the range [0, d). In the literature, there exist (at least) 
two conventions for the gradient and the divergence operators depending on whether the spatial 
derivative axis is appended or prepended; or informally, whether gradients and divergences are 
taken column-wise or row-wise. We denote the former convention by grad(w) and the latter 
by V«. The two choices are reflected in UFL via two gradient operators: grad(A) corresponding 
to grad(u), and nabla_grad(A) corresponding to Vv. The two divergence operators div(A) and 
nabla_div(A) follow the corresponding traditions. Also available are the curl operator and its 
synonym rot, as well as a shorthand notation for the normal derivative. 

Expressions can also be differentiated with respect to user-defined variables. An expression e 
is annotated as a variable via v = variable (e) . Letting A be some expression of v, the derivative 



Math, notation 


UFL notation 


Math, notation 


UFL notation 






A l or |± 

dA 
dn 


A.dx(i) or Dx(A, i) 


(T, ifc 


conditional (c 




T, F) 






Dn(A) 


1 F, otherwise 




divA 


div(A) 


a = b 


eq(a, b) 






grad ^4 


grad(A) 


a^b 


ne(a, b) 






V ■ A 


nabla_div(A) 


a<b 


le(a, b) or a 


<= 


= b 


VA = V <%>A 


nabla_grad(A) 


a>b 


ge(a, b) or a 


>= 


= b 


curl A = V x A 


curl (A) 


a <b 


It (a, b) or a 


< 


b 


rot A 


rot (A) 


a> b 


gt (a, b) or a 


> 


b 


v — e 


v = variable (e) 


lAr 


And(l, r) 






dA 


diff(A, v) 


I V r 


0r(l, r) 






df 


exterior_derivative (f ) 


->c 


Not(c) 







Table 8: (Left) Table of differential operators. (Right) Table of conditional operators. 



12 



Mathematical notation 


UFL notation 


/+,/" 


f ('+'), f ('-') 


(f) 


avg(f ) 


ill 


jump(f) 


Uln 


jump(f, n) 



Mathematical notation UFL notation 



>Tk- 



1,2 



Idx 
Ids 
Ids 



h+l 



I * dx(k) 
I * ds(k) 
I * dS(k) 



Table 9: (Left) Table of discontinuous Galerkin operators. (Right) Table of subdomain integrals. 



of A with respect to v then reads as diff (A, v). Note that diff (e, v) == 0. 

3.2.6 Conditional operators 

The lack of control flow statements and mutable variables in UFL is offset by the inclusion of 
conditional statements, which are equivalent to the ternary operator known from, for example, 
the C programming language. To avoid issues with overloading the meaning of some logical 
operators, named operators are available for all boolean UFL expressions. In particular, the 
equivalence operator == is reserved in Python for comparing if objects are equal, and the delayed 
evaluation behavior of Python operators and, or and not preclude their use as DSL components. 
The available conditional and logical operators are listed in Table[8] (right) and follow the standard 
conventions. 

3.2.7 Discontinuous Galerkin operators 

UFL facilitates compact implementation of discontinuous Galerkin methods by providing a set 
of operators targeting discontinuous fields. Discontinuous Galerkin discretizations typically rely 
on evaluating jumps and averages of (discontinuous) piecewise functions, defined relative to a 
tessellation, on both sides of cell facets. More precisely, let F denote an interior facet shared by 
the cells T + and T~ , and denote the restriction of an expression / to T + and T~ by /+ and /~, 
respectively. In UFL, the corresponding restrictions of an expression f are expressed via f (' + ') 
and f ( ' - ' ) . 

Two typical discontinuous Galerkin operators immediately derive from these restrictions: the 
average (/) = (/+ + f~)/2, and jump operators [/] = /+ — /~. For convenience, these two 
operators are available in UFL via avg(f) and jump(f). Moreover, it is common to use the 
outward unit normal to the interior facet, n, when defining the jump operator such that for a 
scalar-valued expression [/]„ = f + n + + f~n~, while for a vector- or tensor-valued expression 
[/In = / + ' n+ + / _ ' n - These two definitions are implemented in a single UFL operator 
jump(f , n) by letting UFL automatically determine the rank of the expression and return the 
appropriate definition. The available operators are presented in Table [9] (left). 

Since UFL is an embedded language, a user can easily implement custom operators for dis- 
continuous Galerkin methods from the basic restriction building blocks. The reader is referred 
to 01gaard et al. 37J for more details on discontinuous Galerkin methods in the context of a 
variational form language, developed for FFC, and automated code generation. 

3.3 Integrals and variational forms 

In addition to expressions, UFL provides concepts and syntax for defining integrals over facets and 
cells and, via integrals, for defining variational forms. Variational forms can further be manipulated 
and transformed via form operations. This sublanguage is described in the sections below. 

3.3.1 Integrals and forms 

The integral functionality provided by UFL is centered around the integrals featuring in variational 
forms as summarized by the canonical expression ([8]). An integral is generally defined by an 
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Mathematical notation UFL notation Form type 



Ml f-)= L- 1 ~C dx 
±y± \J ' i jji i+/ 2 


M = 


(1-f **2)/(l+f **2) * dx 


Functional 


M(f, g; ) = J Ti grad / • grad g dx 


M = 


dot (grad (f) , grad(g)) * dx 


Functional 


L(v) — Jj- sui(ttx)v dx 


L = 


sin(pi*x [0] ) *v * dx 


Linear form 


L (g; v ) = (g rad g-n)vds 


L = 


Dn(g)*v * ds 


Linear form 


a(u, v) = Jj. uv dx 


a = 


u*v * dx 


Bilinear form 


a(u, v) — Jj- uv dx + fuv ds 


a = 


u*v * dx(0) + f*u*v * ds(l) 


Bilinear form 


a{u,v) — Jj- a (u)(v) ds 


a = 


avg(u) *avg(v) * dS 


Bilinear form 


a(A; u, v) — Jj. AfjU^Vj dx 


a = 


A[i,j]*u.dx(i)*v.dx(j) * dx 


Bilinear form 



Table 10: Table of various forms with reference to coefficient functions f , g; argument functions 
u, v; predefined integration measures dx, ds, and dS; and subdomain notation as introduced in 
Section |2.1.1[ Recall that the predefined UFL measures dx, ds and dS default to dx(0), ds(0) 
and dS(0), respectively, and that the mathematical notation starts counting at 1, while the code 
starts counting at 0. 



integration domain, an integrand and a measure. UFL admits integrands defined in terms of 
terminal expressions and operators as described in the previous section. The mathematical concept 
of an integration domain together with a corresponding measure is however embodied in a single 
abstraction in UFL, namely the Measure class. This abstraction was inherited from the original 
FFC form language. 

A UFL measure (for simplicity just a measure from here onwards) is defined in terms of a 
domain type, a domain identifier and, optionally, additional domain meta-data. The allowed 
domain types include "cell", "exterior_f acet" and " interior _facet". These domain types 
correspond to Lebesgue integration over (a subset of) cells, (a subset of) exterior facets or (a subset 
of) interior facets. The domain identifier must be a non-negative integer (the index k in Q). For 
convenience, three measures are predefined by UFL: dx, ds and dS corresponding to measures over 
cells, exterior facets and interior facets, respectively, each with the default domain identifier 0. 
However, new measures can be created directly, or by calling a measure with an integer k yielding 
a new measure with domain identifier k. Several other, less common, domain types are allowed, 
such as "surface", "point" and "macro_cell", and new measures are easily added. We refer to 
the UFL documentation for the complete description of these. 

A UFL integral is then defined via an expression, acting as the integrand, and a measure object. 
In particular, multiplying an expression with a measure yields an integral. This is illustrated in 
Table [9] (right), with k denoting the domain identifier. We remark that all integrand expressions 
featuring in an integral over interior facets must be restricted for the integral to be admissible. 
The integrand expression will depend on a number of distinct argument functions: the number of 
such is labeled the arity of the integral. 

Finally, a UFL form is defined as the sum of one or more integrals. A form may have integral 
terms of different arities. However, if all terms have the same arity p we term this p the arity of 
the form. Forms are labeled according to their arity: forms of arity are called functionals, forms 
of arity 1 are called linear forms, and forms of arity 2 are called bilinear forms. We emphasize that 
the number of coefficient functions does not affect the arity of an integral or a form. Table [10| shows 
simple examples of functionals, linear forms and bilinear forms, while more complete examples are 
given in Section |4j 
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Mathematical notation 


UFL notation 


Description 




L = lhs(F) 


Extract terms of arity 2 from F 




a = rhs(F) 


Extract terms of arity 1 from F 




(a, L) = system(F) 


Extract terms of arity 2 and 1 from F 


a i y Oj 


a_star = adjoint (a) 


Derive adjoint form of bilinear form a 


F(f;-)^F(g;-) 


G = replace(F, {f:g}) 


Replace coefficient f with g in F 


F(;-)^F(f;-) 


M = action(F, f) 


Replace argument function 1 in F by f 


F(f;-)^D f F(f;-)[v] 


dF = derivative (F, f, v) 


Differentiate F w.r.t f in direction v 



Table 11: Overview of common form operators, with reference to a given form F of possibly mixed 
arity, coefficient functions f , g and an argument function v. Note that all form operators return 
a new form as the result of the operation; the original form is unchanged. 



Mathematical operation UFL notation 



Li(u,p;v) 


d 

(It 


[M{u + 


TV 'P)] T =0 


LI 


= derivative (M, 


u, v) 


L2(u,p; (?) 


d 

dr 


[M{u,p 


+ r=0 


L2 


= derivative (M, 


p. q) 


L 3 (u,p;w) 


d 

dr 


[M(u-f 


TV,p + Tq)] T=Q 


L3 


= derivative (M, 


(u, p) , w) 


L 4 (u,p; s) 


d 

dr 


[M(u + 


rse y ,p)] T=Q 


L4 


= derivative (M, 


u[l] , s) 



Table 12: Example derivative calls on a functional M : W = V X Q — > K. for a vector- valued 
function space V, with a scalar subspace V\, and a scalar- valued function space Q, and with 
reference to coefficient functions u G V, p G Q, and argument functions v G V, q G Q, s G V\ and 
w G W. 



3.3.2 Form operators 

UFL provides a select but powerful set of algorithms that act on forms to produce new forms. An 
overview of such algorithms is presented and exemplified in Table |11| 



The first three operators in Table 11 lhs, rhs and system, extract terms of certain arities 
from a given form. More precisely, lhs returns the form that is the sum of all integrals of arity 2 
in the given form, lhs returns the form that is the negative sum of all integrals of arity 1 and 
system extracts the tuple of both the afore results; that is, system(F) = (lhs(F), rhs(F)). 
These operators are named by and typically used for extracting 'left-hand' and 'right-hand' sides 
of variational equations expressed as F(u,v) — 0. Next, the adjoint operator acts on bilinear 
forms to return the adjoint form a*: a(;u,v) i— >• a*(;u, v) — a(;v,u). The replace operator 
returns a version of a given form in which given functions are replaced by other given functions. 
The action operator can be viewed as special case of replace: the argument function with the 
lowest count is replaced by a given function. 

The Gateaux derivative ^ is provided by the operator derivative, and is possibly the most 



important form operator. Table 12 provides some usage examples of this operator. With reference 
to Table [l2j we observe that forms can be differentiated with respect to coefficients separately 
[Li and L 2 ) or with respect to simultaneous variation of multiple coefficients (L 3 ). Note that 
in the latter case, the result becomes a linear form with an argument function in the mixed 
space. Differentiation with respect to a single component in a vector-valued coefficient is also 
supported (£4). 

The high-level operations on forms provided by UFL can enable the expression of algorithms 
at a higher abstraction level than what is possible or practical with a traditional implementation. 
Some concrete examples using UFL and operations on forms can be found in [151 SI] • 
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4 Examples 



We present now a collection of complete examples illustrating the specification, in UFL, of the 
finite element discretization of a number of partial differential equations. We have chosen prob- 
lems that compactly illustrate particular features of UFL. We stress, however, that UFL is not 
limited to simple equations. On the contrary, the benefits of using UFL are the greatest for 
complicated, non-standard equations that cannot be easily or quickly solved using conventional 
libraries. Computational results produced using UFL have been published in many works by 
third-parties and by the developers of UFL. Examples of the use of UFL for challenging, complex 
problems from a range of fields, and which make available complete source codes, can be found in 
the references [HHSl HO [33 EH W\ ■ 



4.1 Poisson equation 

As a first example, we consider the Poisson equation and its discretization using the standard 
if ^conforming formulation, a L 2 -conforming formulation and a mixed i?(div)/L 2 -conforming 
formulation. The Poisson equation with boundary conditions is given by 

— div(ftgrad it) = / in f2, u — uq on To, ~nd n u = g on Tn, (20) 

where TdUFn = dQ is a partitioning of the boundary of f2 into Dirichlct and Neumann boundaries, 
respectively. 

4-1.1 H 1 -conforming discretization 



The standard H 1 -conforming finite element discretization of (20) reads: find u 6 V/, such that 



a(u 7 v) = / k grad u • grad v dx — / fvdx — gv ds = L(v) (21) 
Jn Jn Jr N 

for all v € Vh, where Vh is a continuous piecewise polynomial trial space incorporating the Dirichlet 
boundary conditions on Trj and V/, is a continuous piecewise polynomial test space with zero trace 
on Trj . We note that as a result of the zero trace of the test function v on Trj , the boundary integral 



in (21) may be expressed as an integral over the entire boundary dQ. A complete specification 



of the variational problem (21) in UFL is included below. The resulting expression trees for the 
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bilinear and linear forms are presented in Figure [T] 

UFL code 

element = FiniteElement (" Lagrange " , triangle, 1) 

u = TrialFunct ion ( element ) 
v = Te st Funct ion ( element ) 

f = Coeff i c ient ( element ) 
g = Coeff i c ient ( element ) 
kappa = Coef f i c ient ( element ) 

a = kappa * inner ( grad (u) , grad (v) ) *dx 
L = f*v*dx - g*v*ds 
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Form a 



Form L 



Cell integral Cell integral Exterior facet integral 



/L R 



kappa 



R \L 



L \R 



grad 



grad 



Figure 1: Expression trees for the H 1 discretiza- 



tion of the Poisson equation ( 21 ). 



4-1-2 L 2 - conforming discretization 

For the L 2 discretization of the Poisson equation, the standard discontinuous Galerkin/interior 
penalty formulation [TOl reads: find u <E Vh such that 



£ /.grad.-grad.d, 



f eJ™ 



re grad u) ■ \v\ n ds - / (re grad v) ■ \u} n ds 
! Jf 

Kd„uvds— / Kd„vuds+ / — — uvds 



7(k) 



'«1M ds 



fv dx 



gv ds 



d n uov ds — d n vuo ds 
r D ^r D 



1 K 



uqv ds 



(22) 



for all v € Vh, where [u], [w]„ and (w) are the standard jump, normal jump and average operators 



(see Section 3.2.7) 



The corresponding implementation in UFL is shown below. The relevant operators for the 
specification of discontinuous Galerkin methods are provided by the operators jumpO and avg(). 
We also show in Figures [2] and [3] the expression tree for the bilinear form of the L 2 -discretization, 
which is notably more complex than for the -H^-discretization. 



UFL code 
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element = FiniteElement (" Discontinuous Lagrange", triangle, 1) 

u = TrialFunct i on ( element ) 
v = Te st Funct ion ( element ) 

f = Coeff i c ient ( element ) 
g = Coeff i c ient ( element ) 
kappa = Coef f i c ient ( element ) 
uO = Coef f i cient ( element ) 

h = 2 * tr i angle . c ir cumr adius 
n = triangle . n 

gamma = 4 

a = kappa * dot ( grad (u) , grad(v))*dx \ 

- dot ( avg ( kappa * grad (u )) , jump(v, n))*dS \ 
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dot ( avg ( kappa * gr ad (v )) , jump(u, n))*dS \ 
( gamma * avg (kappa ) / avg (h) ) * j ump (u ) * j ump ( v) * dS \ 
dot ( kappa* gr ad (u) , v*n)*ds(0) \ 



- aot uappa*graou; , v*n;*as 

- dot ( kappa * gr ad (v) , u*n)*ds 
+ ( gamma * kappa/h ) *u*v* ds ( ) 



(0) \ 



L = f*v*dx - g*v*ds(l) \ 

- dot ( kappa* gr ad ( uO ) , v*n)*ds(0) - dot (kappa*grad (v) , u0*n)*ds(0) \ 
+ (gamma*kappa/h) *u0*v*ds (0) 



4-.1.3 H (div) /L 2 - conforming discretization 

Finally, we consider the discretization of the Poisson equation with a mixed formulation, where 



the second-order PDE in ( 20 ) is replaced by a system of first-order equations: 



diver = / in tt, cr + Kgradu = in f2, u — uq on Td, a ■ n = g on Tn- (23) 



A direct discretization of the system ([23]) with, say, continuous piecewise linear elements is un- 
stable. Instead a suitable pair of finite element spaces must be used for a and u. An example of 
such a pair of elements is the i?(div)-conforming BDM (Brezzi-Douglas-Marini) element [15] of 
the first degree for a and discontinuous piecewise constants for u. Multiplying the two differential 



equations in (23 1 with test functions r and v, respectively, integrating by parts and summing the 
variational equations, we obtain the following variational problem: find (cr, u) 6 Vh = V£ x V^ 
such that 

/ (div cr)?; dx + / a -rdx— / udiv(nT)dx~ / fvdx — ku t ■ nds (24) 
Jn Jn Jn Jn Jv D 

for all (t, v) G Vh- Note that the Dirichlet condition u = uq on Td is a natural boundary condition 
in the mixed formulation; that is, it is naturally imposed as a weak boundary condition as part of 
the variational problem. On the other hand, the Neumann boundary condition a ■ n = g must be 
imposed as an essential boundary condition as part of the solution space Vh- The corresponding 
specification in UFL is given below. 

UFL code 



9 
10 
11 
12 
13 
14 
15 
16 



BDM = FlnlteElement ("BDM" , triangle, 1) 
DG = FiniteElement ("DG" , triangle, 0) 
CG = FiniteElement (" CG " , triangle, 1) 
element = BDM * DG 

(sigma, u) = Tr i alFunct ions ( element ) 
(tau, v) = Te st Funct i ons ( element ) 

f = Coefficient (CG) 
kappa = Coef f i c ient ( CG ) 
uO = Coef f icient (DG) 

n = triangle . n 

a = div ( sigma) *v*dx + dot(sigma, tau)*dx - u* div ( kappa* tau ) *dx 
L = f*v*dx - kappa* uO * dot ( tau , n)*ds 



4.2 Stokes equations 

As a second example of a mixed problem, we consider the Stokes equations given by 

— Au + gradp = / in fi, divu = in ft, (25) 

together with a suitable set of boundary conditions. The first equation is multiplied with a test 
function v, the second equation is multiplied with a test function q, and after integration by parts, 
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Figure 3: Detail of the 
expression tree of Figure [2] 
for the L 2 discretization of 
the Poisson equation ( [22] ) . 
The expression tree rep- 
resents the expression 
gamma*avg (kappa) / avg (h) . 



the resulting mixed variational problems is: find (u,p) € Vh such that 

/ grad u ■ grad v dx — / (divv)pdx + / (div u)q dx — / fvdx. 
Jn Jn Jn Jn 



(26) 



The corresponding specification in UFL using an inf-sup stable \P2\ d ~P\ discretization (Taylor- 
Hood element [45J) is shown below. 



UFL code 



1 


P2 = VectorElement (" Lagrange " , 


triangle 


, 2) 


2 


PI = FiniteElement (" Lagrange " , 


triangle 


, 1) 


3 


TH = P2 * PI 






4 








5 
6 


(u, p) = TrialFunctions (TH) 
(v, q) = TestFunctions (TH) 






7 
8 


f = Coefficient (P2) 






9 








10 


a = inner ( grad (u) , grad (v) ) *dx 


- div ( v) 


*p*dx + div(u)*q*dx 


11 


L = dot(f, v)*dx 







4.3 Neo-Hookean hyperelastic model 

We consider next a hyperelastic problem posed in terms of the minimization of potential energy. 
For a body f2 C M. d in a reference configuration, where 1 < d < 3, the total potential energy II 
reads 



Jn 



(27) 



B ■ v dx — / T • v ds, 
>n Jr 

where tp is the stored energy density, v is the displacement field, B is the nominal body force and 
T is the nominal traction on the domain boundary T = dfl. For the compressible neo-Hookcan 
model, the strain energy density reads 



X 



V>(«) = f(/c-3)- M lnJ+-(lnJ) 2 , 



(28) 



where J is the determinant of the deformation gradient F = I + grad u and I c is the trace of the 
right Cauchy-Green tensor C — F T F. Solutions u to the hyperelastic problem minimize (27): 



u — argmin II(u) 



(29) 
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The classical variational approach to solving (29 1 involves taking the first variation/Gateaux 
derivative of (27 1 with respect to v and setting this equal to zero for all v, which in turn can 



be solved using Newton's method by taking a second Gateaux derivative to yield the Jacobian. 
The corresponding specification in UFL for such an approach is shown below. 

UFL code 



# Function spaces 

element = VectorElement (" Lagrange " , tetrahedron, 1) 



# Trial and test functions 
du = TrialFunc t i on ( element ) 
v = Te st Func t i on ( e lement ) 

# Functions 

u = Coeff i c ient ( element ) 
B = Coeff i c ient ( element ) 
T = Coeff icient ( element ) 



# Incremental displacement 

# Test function 



# Displacement from previous iteration 

# Body force per unit volume 

# Traction force on the boundary 



# Kinematics 

I = Identity ( element . cell (). d) # Identity tensor 
F = I + grad(u) 
C = F.T*F 



# Deformation gradient 

# flight Cauchy-Green tensor 



# Invariants of deformation tensors 
Ic , J = tr (C) , det (F) 

# Elasticity parameters 

mu = Const ant ( t etrahedr on ) 

lmbda = Const ant ( t etrahedr on ) 

# Stored strain energy density (compressible neo -Hookean model) 
psi = (mu/2)*(Ic - 3) - mu*ln(J) + ( lmbda/2 ) * ( In ( J ) ) * * 2 

# Total potential energy 

Pi = psi*dx - inner(B, u)*dx - inner(T, u)*ds 

# First variation of Pi (directional derivative about u in the direction of v) 
F = derivative (Pi , u, v) 

# Compute Jacobian of F 
J = der i vat i ve (F , u, du) 



A complete solver using precisely the above formulation is distributed as a demo program as 
part of DOLFIN 0H1]. 



4.4 Constrained optimization 

As a final example, we consider a PDE constrained optimization problem. We wish to determine 
the control variable p that minimizes the cost functional 

J(u,p) = [ -(u — u) 2 dx + I -ap 2 dx, (30) 
Jn 2 Jn 2 

where u is a given target function, a is a regularization parameter and u is the state, constrained 
by the variational problem 

a(u,v) = / uv -\- grad u • grad v dx — / pv dx = b(p; v), (31) 
Jn Jn 

which should hold for all test functions v in some suitable test space. 

A way to solve the optimization problem is to define the Lagrangian functional C as the sum 
of the cost functional and the weak constraint, in which v plays the role of a Lagrange multiplier: 

C(u,p, v) = J{u,p) + a(u,v) - b(p; v). (32) 
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The solution of the optimization problem can be found by seeking the stationary points of the 
Lagrangian: 

D u C[u,p,v)[8u]= V 6u, (33) 
D p C(u,p,v)[5p}= V 8p, (34) 
D v C{u,p,v)[5v}=0 y Sv, (35) 

which in this case is a linear system of equations. The corresponding implementation in UFL is 
shown below. This example shows how to export a specific set of forms from a .ufl file. Without 
the explicit forms = [. . .] line, UFL will export forms with predefined names, namely a, L, F, J 
and M. 



UFL code 



1 


# Define mixed function space for all variables 


2 


V = FiniteElement (" Lagrange " , triangle, 1) 


3 
4 


W = MixedElement (V , V, V) 




5 


# Define coefficients 




D 


w = Coefficient (W) # 


mixed function with all unknowns 


7 


u = w[0] # 


state 


8 


p = w[l] # 


control 


9 


v = w[2] # 


Lagrange multiplier 


10 


alpha = Coef f icient (V) # 


regular ization parameter 


11 


ubar = Coef f i c i ent ( V) # 


observation 


12 


pbar = Coef f i c i ent (V) # 


analytic control 


13 






14 


# Define forms 




15 


def J(u, p) : return 0.5*(u 


- ubar)**2*dx + . 5 * alpha *p* * 2 * dx 


16 


def a(u, v) : return (u*v + 


dot(grad(u), grad(v)))*dx 


17 


def b(p, v) : return p*v*dx 




18 






19 


# Define Lagrangian 




20 


L = J(u, p) + a(u, v) - b(p, v) 
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22 


# Differentiate 




23 


F = derivative (L , w) 




24 


dF = derivative (F , w) 




25 


mF = -F 




26 






27 


# Define some norms 




28 


L2p = . 5* (p-pbar ) **2*dx 




29 


L2u = . 5* (u-ubar ) **2*dx 




30 


J = J(u, p) 




31 






32 


# Export forms 




33 


forms = [mF , dF , J, L2p , L2u] 



5 Representation of expressions 

UFL is a collection of value types and operators. These types and operators have been presented 
in the preceding sections from a user perspective, with a focus on their mathematical definitions. 
To discuss the representation and algorithms in the symbolic framework underlying the language 
implementation, we will now pursue a more abstract approach. 

5.1 Abstract categorization of expression types and DAG 

As a domain specific language embedded in Python, UFL does not have a formal grammar of its 
own. However, it is still useful to compactly categorize the various types of expressions before 
proceeding. An expression e can be either a terminal expression t, or an operator o applied to one 
or more other expressions. A terminal expression t can be categorized in one of the groups: 
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i Multi-index. 



I Literal constant. 

g Geometric quantity. 

/ Function. Further classified as: 

c Coefficient function. 
a Argument function. 

An operator type o can be categorized in one of the following groups: 

/ Indexing related operators, which manipulate the indices and shape of expressions. These 
operations do not introduce any computation. 

A Basic arithmetic operators, including the operators +, -, * and /. 

F Nonlinear scalar real functions, such as ln(x), exp(a;), sin(x) and arccos(x). 

T Tensor algebraic operators. These are convenience operators acting on nonscalar expressions, 
such as the dot, inner, and outer products, and the transpose, determinant or deviatoric part. 

D Differential operators, including both spatial derivatives, derivatives with respect to expres- 
sions annotated as variables and directional derivatives with respect to coefficient functions. 

R Restrictions of functions to cells and related operators such as the jump and average between 
cells. 

B Boolean operators, including =, 7^, <, >, <, >, A, V, -1. These can only be used within a 
conditional operator. 

C The conditional operator: the ternary operation "/ if condition else g" . This is the only 
explicit branching instruction, not counting restrictions. 

With these definitions, the following diagram gives a compact semi-formal overview of the types 
and operators that an expression e can be recursively built from: 



As before, t represents any terminal expression, while o represents any operator. If an expression 
e is an application of an operator, we say that it depends on the expressions e\, . . . , e n . Depen- 
dencies between expressions are, by construction, one-way relations and can therefore be viewed 
as the edges of a directed acyclic graph (DAG) , where each vertex is itself an expression e, . Each 
graph vertex has an associated type, whether it is a terminal type or operator type. Vertices rep- 
resenting a terminal expression may have additional descriptive data associated with them, while 
vertices representing non-terminal subexpressions are fully identified by the type and sequence of 
child vertices. We note here that the ordering of child vertices (dependencies) is important (for 
nonsymmetric operators) . 



o 



e 



t 



[i\l\g\c\a] 

[I\A\F\T\D\R\B\C] 
[i|o(ei, . . . ,e n )\ 
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5.2 Representation of directed acyclic graphs 

From an implementation viewpoint, a DAG can be represented by link-based or list-based data 
structures. In a link-based data structure, each vertex of the graph is represented by a single typed 
object. Objects of operator types store references to the objects representing their operands. These 
references are the edges of the DAG. Objects of terminal types may store additional data. This is 
the natural DAG representation for a symbolic library while an expression is being built, and is 
the way UFL expressions are represented most of the time. When a UFL operator is evaluated in 
Python, the result is a new object)^] representing the operator. Thus, both the time and storage 
cost of building an expression is 0(n) in the number of operators applied n. Since each operator 
invocation leads to a new graph vertex without any global knowledge available, some duplicate 
expressions may occur in this DAG representation. 

A list-based data structure for an expression DAG can be constructed from the link-based data 
structure, if needed. In this data structure, references to each unique subexpression are placed in 
a topologically sorted vertex list and all edges (dependencies) are stored in another list as a tuple 
of vertex list indices. While the vertex list is constructed, duplicate expressions with the exact 
same symbolic representation are automatically mapped to the same vertex. This mapping can 
be achieved through 0(1) insertion into a hash map, which retains the overall 0(n) performance 
of building this data structure. The main advantages of the list-based graph representation are 
efficient ordered iteration over all vertices and easy access to dependency information. This is 
beneficial for some algorithms. 

5.3 Expression type class hierarchies 

Each DAG vertex, or expression object, is an instance of a subclass of Expr. The type of terminal 
or operator is determined by the subclass, from a class hierarchy which is divided into Terminal 
and Operator subtypes. Concrete Operator subclasses store references to the objects representing 
their operands. Since a non-terminal expression is uniquely determined by its type and operands, 
any other data stored by these classes is purely for performance or convenience reasons. Concrete 
Terminal subclasses, however, can take any necessary auxiliary data in their constructor. This 
data must still be immutable, such that any expression object is hashable. Most operator over- 
loading is applied directly to the Expr class. For example, adding any two expressions will result 
in a new instance of a Sum. All expression classes must overload a set of basic methods. One 
such method is operands (), which returns a tuple of expression objects for the operands of an 
Operator, or the empty tuple for a Terminal object. We refer to the UFL source code [4] for 
further internal implementation details. 

5.4 Representation of indexed expressions 

In some symbolic libraries, expressions with shapes or indices are built on top of an otherwise 
mainly scalar framework. In contrast, UFL considers shapes and indices as an integral part of 
any expression. To support this, expressions provide three methods, namely, shape () which re- 
turns a tuple of positive integers representing the value shape of the expression, f ree_indices () 
which returns an unordered tuple of Index objects for each free-index in the expression and 
index_dimensions () , which returns a mapping from Index objects to the dimensions they im- 
plicitly range over. Together, these methods provide a rich and flexible description of the value 
shape and free-index set. These properties generalize to arbitrarily shaped expressions with arbi- 
trary sets of free-indices, and are accessible for any expression regardless of type. 

In traditional programming languages, the indexing operator [] usually extracts a value from 
a container. In a symbolic environment, the result may instead be an object representing the 
operation of extracting a component. When indexing with free-indices, extraction is not possible 
since the index represents a range of values. When indexing a tensor-valued expression in UFL, 

4 Sometimes more than one new object, but the number of auxiliary objects is always bounded by a small 
constant. 
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the result is represented as an object of the Indexed class, with the original expression and 
a Multilndex object as operands. A Multilndex object^] represents a sequence of Index and 
Fixedlndex objects. The opposite operation from indexing, mapping an expression with free- 
indices to be seen as a tensor-valued expression, is represented with the ComponentTensor class, 
which similarly has the original expression as well as a Multilndex as operands. Expressions 
which contain layers of mapping back-and-forth between index and tensor notation may appear 
complex, but a form compiler can reduce this complexity during the translation to low-level code. 

The implicit summation notation in UFL is applied early in the application of the * operator 
and the derivative operator .dx(i). Subsequently the summation is represented as explicit sums 
over free-indices. This way, the implicit summation rules need no special consideration in most 
algorithms. For example, the sum u^i = UiVi is expressed in UFL code as u[i]*v[i], but is 
represented in the DAG as IndexSum (Product (u[i] , v[i]), (i,))|^] A Product object can be 
constructed directly in algorithms where implicit summation is not the wanted behavior. When 
interpreting an IndexSum object in symbolic algorithms, the range of the sum is defined by the 
underlying summand expression; this is possible as any expression in UFL knows its own shape, 
free-indices and dimensions of its free-indices. 

5.5 Simplification of expressions 

Automatic simplification of expressions is a central design issue for any symbolic framework. At 
one end of the design spectrum are conservative frameworks that preserve expressions exactly as 
written by the user, analogous to an abstract syntax tree of the program input to a compiler. At 
the other end of the spectrum are frameworks where all expressions are transformed to a canonical 
form at construction time. The design of UFL is guided by the intention to be the front end to 
multiple form compilers, and is therefore fairly conservative. If a form compiler does no rewriting, 
the generated code should transparently resemble the UFL expressions as authored by the end- 
user. However, UFL is also a collection of symbolic algorithms, and the performance and memory 
footprint of such algorithms can be significantly improved by certain automatic simplifications of 
expressions at construction time, if such simplifications reduce symbolic expression growth. 

Because form compilers may employ different expression rewriting strategies, we wish to avoid 
doing any simplifications that might remove information a form compiler could make use of. A 
canonical sorting of two sum operands eases the detection of common subexpressions, such as in 
the expression rewriting (b + a)(a + b) — > (a + b) (a + b). However, a canonical sorting of more than 
two sum operands may hide common subexpressions, such as in the rewriting ((a + c) + b)(a + c) — > 
(a + b + c)(a + c) . Therefore, sums and products are represented in UFL as binary operators (as 
opposed to storing lists of terms or factors), with their two operands sorted canonically. Even 
more importantly, we avoid any symbolic rewriting that can lead to numerically unstable floating 
point expressions in generated code. An example of such unsafe operations is the expansion of a 
factored polynomial (a — b)(a — b) — > a 2 — 2ab + b 2 , which becomes numerically unstable in inexact 
floating point arithmetic. 

The performance of the symbolic algorithms and the form compilation processes poses a final 
limitation on the type of automatic simplifications we may apply. It is crucial that the overall 
form compilation process can be designed to have an asymptotic cost of 0(n) in time and memory 
usag^J with n a measure of the length of the integrand expression. Because simplification of 
expressions at construction time is performed once for each expression object, the cost of applying 
the simplification must be a local operation independent of the size of the operand expressions. 
Thus, automatic simplifications which would involve traversing the entire subexpression DAG for 
analysis or rewriting are never attempted. 

5 The type system in UFL is fairly simple, and Multilndex is one of the few exceptions in the Expr type hierarchy 
that does not represent a tensor-valued expression. 

6 Here, the representation of u[i] , v[i] , (i,) are simplified for compact presentation. 

7 Form compilers are free to apply more expensive strategies, but UFL must not render efficient algorithms 
impossible. 
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The following small set of safe and local simplifications are applied consistently when con- 
structing expressions: 



Multiply by one: Ix — > x. A simplification which keeps one operand intact and throws away 
the other one is safe. 

Add zero: + x — > x. 

Multiply by zero: Ox — > 0. To avoid losing the tensor properties of x, we annotate the rep- 
resentation of with the same properties. Therefore we have a special zero representation 
with shape and free-indices. 

Constant folding: f(li,h) —> h- Here h, 1% and I3 are literal constants and / is a function 
or operator that can be computed numerically by UFL. This happens recursively and so 
constant scalar expressions are effectively folded to a single literal constant. 

Canceling indexing: A a E a — > A. The mappings between tensor-valued and indexed expression 
representations cancel when the same multi-index a is used in both operations. If the inner 
and outer multi-indices are not equal, this canceling operation would require rewriting the 
representation of A, which is not a local operation, and therefore not invoked. 

Note that some of these operations will occur frequently during automatic differentiation. For 
example, consider the derivative d(fg)/df, where / and g are independent functions: 

d{fg) df dg 

^T = df g + f W =9 + f0 9 + g - (36) 

Similarly, during some symbolic algorithms, tensor-valued subexpressions are indexed to simplify 
computation and later mapped back to tensor-valued expressions. This process leads to superfluous 
indexing patterns which causes the DAG to grow needlessly. Applying the indexing cancellations 
helps in avoiding this DAG growth. 



6 Algorithms 

The UFL implementation contains a collection of basic algorithms for analysis and transformation 
of expressions. These algorithms include optimized Python generators for easy iteration over 
expression nodes in pre- or post-ordering and iteration over terminal expressions or, more generally, 
expressions of a particular type. In the following, some of the core algorithms and building blocks 
for algorithms are explained. Particular emphasis is placed on the differentiation algorithm. To 
avoid verbose technical details in the discussions of symbolic algorithms, we first define some 
mathematical notation to support abstract algorithm descriptions. For implementation details, 
we refer to the UFL source code [3]. 



6.1 Evaluation algorithm 

Assume that an expression e is represented by a DAG with m terminal and n non-terminal subex- 
pressions, recalling that each subexpression is a DAG vertex. Denote the terminal subexpressions 
by ei, for i = 1, 2, . . . , m and the non-terminal subexpressions by e^, for i = m + 1, to + 2, . . . , m + n. 
For each e^, let P = (^j)jLi be a sequence of pi integer labels referring to the operand expressions 
of et. Note that for each terminal expression, this sequence is empty. Moreover, we require the 
labels to fulfill a topological ordering such that 

Ij<i, j = l,...,pi, i = l,...,m + n. (37) 

We thus have k < i whenever depends on e& (directly or indirectly). Equivalently e; is inde- 
pendent of any subexpression for k > i. 
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We can now formulate two versions of an algorithm for evaluating e = e m+n . Note that these 
algorithms are merely abstract tools for describing the mathematical structure of algorithms that 
follow. For each specific algorithm, we require an evaluation operator E where Vi = E(e.i) is 
called the value of e,. This value can be a floating point value, a new symbolic expression or a 
generated source code string, depending on the purpose of the algorithm. The implementation of 
the evaluation operator E will in general depend on the UFL expression type of its argument e, 
which we denote type(e). 

First, partition the set T of all UFL expression types into disjoint sets of terminal types 
and operator types To, and let Et and Eo be non-recursive evaluation operators for terminal and 
non-terminal expressions, respectively. We can then design a simple recursive evaluation algorithm 
Er of the form 



_ E T {ei), if typefo) G T T) 

£ "" ;J -<- (e.AEuie^u), iftype( ei )G7b. [ > 



E, 



o 



The algorithmic structure of (38) assumes that any subexpression can be evaluated given the 
values of its operands, which is not always true. We therefore partition To further into algorithm- 
specific disjoint sets Ty and Tc, where Ty includes types of expressions that can be evaluated 



given the values of its operands as in ( 38 ) , and Tc includes types of expressions which provide a 



context for the evaluation of the operands. By defining corresponding evaluation operators Ey 



and E c , we can then extend (|38|) into 
E R (ei) 



E T (ei), if type(e 4 ) G T T , 

EyL^ERie^f^y if type(ei) G Ty, (39) 



Ec (e,:, (ej<)?=i > if typefe) G T< 



c- 



Note that the evaluation operator Ec may make further recursive calls to Er or a related recursive 
algorithm, but it is assumed that its operands cannot be pre-evaluated without the provided 
context. The difference between how Ey and Ec are applied corresponds to a post-order versus 
pre-order evaluation of the DAG vertices. 



The recursive operator Er defined by ( 38 1 can be implemented by traversing the link-based 



DAG representation with a post-order traversal algorithm, to recursively visit and evaluate child 



vertices before their parent. The more flexible operator Er defined by (391 can be implemented 
similarly but with a mix of post-order and pre-order traversals depending on the visited types. 
The evaluation of e can then be written simply 

v = E R (e) = E R (e m+n ). (40) 



The simpler recursive operator ( 38 1 can also be implemented as a loop over subexpressions in a 
topological ordering, as shown in Algorithm [T] This non-recursive algorithm can be implemented 
by first constructing the list-based DAG representation and then iterating over the vertices. Using 
the list-based representation and the non-recursive algorithm has the advantage of never visiting 
a vertex twice, even if it is reachable through multiple paths. We remark that a cache mechanism 
may, however, remove such duplicate evaluation in the recursive implementation as well. The 
list-based representation also assigns an integer label to each vertex which can be used to store 
associated data efficiently in arrays during the algorithm. However, the construction of the list- 
based DAG representation is not free, and the fixed labeling is a downside in algorithms where new 
expressions are constructed. In the following exposition, the choice of algorithm structure, recur- 
sive or non-recursive, is considered mainly an implementation detail, controlled by performance 
and convenience considerations and not by functionality. 

6.2 Type based function dispatch and the visitor pattern 

In an implementation of the evaluation algorithm described in the previous section, the specific 
evaluation actions must be selected dynamically based on the type of the expression argument. 
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Algorithm 1 Non-recursive algorithm for evaluation of expressions, 
l: for i = 1, . . . , m do 

2: Vi = E T (ei) 

3: end for 

4: for i = m + 1, . . . ,m + n do 

5: Vi = EoiCi^Vji)^) 

6: end for 

7: U := U m+ „ 



By subclassing the UFL provided class Mult iFunct ion and implementing a handler function for 
each expression type, calls to an instance of this class are dynamically dispatched to the correct 
handler based on the type of the first argument. If a handler is missing, the closest superclass 
handler is used instead, which makes it easy to implement default rules for groups of types. An 
example is shown below. 

Python code 

1 from ufl import triangle 

2 from ufl . algorithms import Mult iFunct i on 

3 class Exampl e Algorithm ( Mult iFunct i on ) : 

4 def t erminal ( self , e): return "unhandled terminal type" 

5 def zero (self , e): return repr(e) 

6 E = Example Algorithm ( ) 

7 print E ( triangle . n ) # Prints "unhandled terminal type" 

8 print E (0* triangle . n) # Prints "Zero ((2 ,) , () , {})" 



Building on this same dynamic multifunction design and the Visitor pattern |17j . the class 
Transformer can be subclassed in the same way to implement many recursive symbolic algorithms 



following the structure of (38) or (39). Calling upon an object of a Transformer subclass to 
visit an expression will result in a recursive application of type-specific rules to subexpressions. 
The example below shows a numerical evaluation of a simple expression, using a pure post-order 



implementation as in ( 38 1 . Whether to visit expressions post- or pre-order is specified per handler 
simply by taking visited expressions for child nodes as arguments or not. In more detail, and with 
reference to the example below, a type will be placed in the Tq set if the corresponding handler 
omits the * values argument. The visit method will then automatically call the handler without 
first handling the operands of its argument. 

Python code 



from ufl import triangle 

from ufl . algorithms import Transformer 

class Example Algorithm ( Transf ormer ) : 

def cell_volume ( self , e): return 0.5 

def scalar .value ( self , e) : return float(e) 

def product ( self , e, *values): return values [0] * values [1] 

def division ( self , e, *values): return values [0] / values [1] 

E = Example Algorithm ( ) 

print E. visit (3 * triangle . volume / 2) # Prints 0.75 



6.3 Partial evaluation 

Some symbolic algorithms involve modification of subexpressions, and such algorithms share a need 
to apply an operator to a new sequence of operands. We will designate the notation type(e)((/,)j) 
to the construction of an operator of the same type as e with the given operands. If the operands are 
unchanged, the original expression can be reused since all expressions are considered immutable, 
thus saving memory. This can easily be accomplished for UFL-bascd algorithms by subclassing 
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the ReuseTransf ormer class. In this case, the algorithm inherits the fallback rules given by 



, if type(ej) e T T , 

L, " : ' { "' i = \tyve(e z ) ((E^iej^u) , if typefc) e T Q , 

where -E S eif refers to -Excuse or overridden rules in a subclass. Based on -Grouse, algorithms can be 



written to just modify what they need and let the fallback rules in (41) rebuild the surrounding 
expression with no additional algorithm-specific code. This allows a very compact implementation 
for algorithms such as the partial evaluation in which terminal expressions are replaced with other 
expressions through a given mapping. As an example, consider a partial evaluation algorithm 
mathematically described by: 

fmap(e f ), if type( ei ) € T T , 

^replace (ti) = < . (42) 

I delegate to-E/ reuse (ei), otherwise, 

where delegate to E reuse (ei) represents delegation to inherited rules from the superclass. This can 
be compactly implemented in UFL as follows: 

Python code 



from ufl . algorithms import ReuseTransf ormer 

class Replacer ( ReuseTransf ormer ) : 

def __ini t __ ( self , mapping): 

ReuseTransf ormer . __init__ (self) 

self . .mapping = mapping # A Python diet 

def terminal ( self , e): 

return self . .mapping . get (e , e) 



6.4 Differentiation 

The differentiation algorithm in UFL is a two-level algorithm. The two levels are used to handle 
expressions involving derivatives with respect to different types of variables. In a first step, a sim- 
ple outer algorithm is employed to evaluate the innermost derivatives first; that is, the derivative 
expressions closest to the terminal expressions. This outer algorithm then calls a single- variable 
differentiation algorithm for each derivative expression visited. This, in turn, allows the inner 
algorithm to assume that no nested derivatives are encountered. Thus, for instance in the eval- 
uation of d(cgrad(vu))/dv, an inner algorithm is called first to evaluate grad(i>it), and second to 
evaluate d(c(grad(w)u + vgrad(u)))/dv. Note that the grad operator is kept in the expression 
DAG after derivative evaluation, but is guaranteed to only apply directly to spatially varying 
terminal expressions. More sophisticated approaches to nested differentiation have been explored 
in Karczmarczuk [19] , Pearlmutter and Siskind [38] and Siskind and Pearlmutter [44] , however we 
have considered this additional complexity unnecessary for the purpose of UFL. 

The inner algorithm handles differentiation of an expression e with respect to a single differ- 
entiation variable u. Both e and u may be tensor-valued, but in the following we illustrate using 
scalars for the sake of clarity. By setting the value in Algorithm [T] to a tuple of the subexpression 
and its derivative: Uj = (ei,dei/du), we obtain the standard Forward-mode Automatic Differ- 
entiation algorithm^ (see [18]). This algorithm can also be written in recursive form as in (38) 



However, because of a few exceptions specific to differentiation variable types, the actual algorithm 



in UFL requires the more flexible framework given by ( 39 ) . For simplicity, we will use = de^ / du 
to define the evaluation rules below. 



8 However, output of the algorithm is a new symbolic UFL expression, and the algorithm is therefore clearly a 
symbolic differentiation algorithm. In UFL context, the distinction between Automatic Differentiation and Symbolic 
Differentiation is therefore irrelevant. 
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Generic differentiation rules arc implemented as handler functions in a transformer class cor- 
responding to 



0, if type(ei) G T T , 

^AD(ei) - { e p Li a^ EAn{eij)j a type(ei) e Tq 



(43) 



For each type of differentiation variable, the default differentiation rules in -Ead are subclassed to 
encode the dependency of expression types with respect to the differentiation variable. 

For spatial derivatives, the full gradient is used to represent the derivatives of functions, giving 
the evaluation operator: 



I, 
0, 

grad e, 
grad grad /, 



if e is the spatial coordinate vector, 
if e is a piecewise constant function, 
if e is a non-constant function, 
if e is grad /, 



(44) 



delegate to Ead (e) , otherwise. 



For differentiation with respect to user-annotated variables (see 3.2.51, the operator rules are 
instead modified as: 



1, if e is the variable instance?/, 

-EVd(c) = \ Evd(/) 7 if e is another variable instance annotating the expression /, (45) 

delegate to -Ead (e) , otherwise. 

Note that in this case the Variable type lies in the Tc set, and is thus visited before the annotated 
expression. Hence, the underlying expression is visited by a further recursive call to Evd(/) only 
if the variable e is different from u. 

Finally, for directional derivatives of an expression e with respect to a coefficient function u in 
the direction of v, with possibly user-specified dg/du = h for a subexpression g, the rules become: 



E DD (e) = < 



v, if e is the function u, 

grad E DD (/), if e is grad /, 

hv, if e is the function;;, 

delegate to Ead (e), otherwise. 



(46) 



Note that in this case, we consider the gradient of a function as a terminal entity, although it is 
represented as two expression nodes. 

To support differentiation with respect to specific components of a mixed function (or a vector- 
valued function), the same rules can be applied by choosing an appropriate expression for v. As 
an example, consider the functions 



u : X 



u = (ui,u 3 ), 



v : X 



(47) 



To compute D^M{u)[v\, the differentiation rule for u must yield a vector- valued derivative of 
the same value shape as u. This is accomplished by padding v and using v — (#1,0,1)2) as the 
direction. That is, 



DuM (u) [v] = -j- [M ((wi + TV! , u 2 + Or, u 3 + tv 2 ))] t=q = D U M (u) [v] 



(48) 



This concept of padding to support component-wise derivatives extends to arbitrary tensor- valued 
functions and functions in mixed element spaces, as well as differentiation of variable components 
that are part of different functions. 
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We conclude this section by commenting on the relation between the UFL differentiation 
algorithms and the algorithms implemented in Sundance |28j . In the Sundance approach, the 
derivative is computed numerically on the ffy, and the use of BLAS amortizes the cost of expression 
DAG traversal at run time. However, the traversal cost does increase with the size of the expression 
DAG. Therefore, Sundance avoids computing the DAG for the derivative. In contrast, as UFL is 
typically used in combination with code generation tools, we may differentiate and then simplify. 
This allows us to produce the expression DAG for the derivative and then produce efficient code 
for it. We can see the symbolic traversal cost as a part of the software build time, and it does not 
affect the runtime for computing/assembling variational forms. 



7 Validation 



Several validation steps are performed by UFL at the stage where an operator is applied to 
an expression. All UFL operator types validate the properties of its operands in various ways to 
ensure that the expressions are meaningful. Most importantly, each operator validates the operand 
value shapes and verifies that the use of free-indices is consistent. This type of validation catches 
common indexing bugs at an early stage. Other examples of validations include: checking that 
value shapes and indices match when adding expressions; checking value shapes for tensor algebra 
products; and checking of index ranges for explicit indexing of tensor-valued expressions. Most 
indexing in UFL expressions use implicit ranges, which reduces repetition and common sources of 
errors. When defining a form, an integrand expression must always be scalar-valued without any 
free-indices. This is also checked. When the form is preprocessed, prior to form compilation, a 
number of properties are checked such as: all integrals must depend linearly on the same set of 
argument functions; and in interior facet integrals all functions must be restricted. The latter in 
particular forces increased clarity in the formulations. 

Testing of symbolic frameworks is hard because every algorithm must be tested with an ap- 
propriate selection of expression type combinations to achieve high test coverage. In an attempt 
to answer to this challenge in UFL, we have used multiple layers of defensive programming with 
assertions, unit testing and integration testing in other FEniCS components. Static code analysis 
with PyChecker [3] was very useful during the main development phase. Algorithms in UFL are 
sprinkled with assertions to document assumptions and catch any that fail. Unit tests cover many 
common (and uncommon) combinations of operators and application of algorithms such as differ- 
entiation. In addition to the unit tests in UFL, unit and regression tests in FFC and DOLFIN 
test the use of UFL for many PDE application examples. 

As an example for testing differentiation, a numerical evaluation of the symbolic derivative can 
be compared with the evaluation of manually derived derivative expressions or reference values. 
Of particular interest when considering validation of UFL is a set of integration tests in DOLFIN 
which exploit Green's theorem in ID or nD for a scalar function f(x) or a vector-valued function 
v(x), that is 

f'(x) dx — [f(x)] b , divvdx = v-nds. (49) 
Jn Jdn 

This identity is ideal for combined testing of symbolic differentiation, numerical differentiation, 
and symbolic evaluation, or even higher order derivatives by setting v = grad / in ( 49 1 . Exploiting 
such mathematical identities is key to robust testing of mathematical software, and combining the 
symbolic and numerical paradigms provides good opportunities for discovering errors. 

A software stack such as that provided in the FEniCS Project with (i) a DSL and symbolic 
framework in UFL, (ii) automated code generation in the form compiler FFC [26], and (iii) library 
code in the problem solving environment DOLFIN |27| , naturally introduces additional complexity 
to the debugging process if something goes wrong. However, the high abstraction level allows us 
to check the correctness of end-user programs in various ways. Automated validation of expres- 
sions and forms by UFL allows consistency checks and catching of user errors at various levels of 
abstraction. Last, but not least, UFL has been tested in active use by researchers for more than 
three years as part of the FEniCS Project. 
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8 Conclusions 



We have presented the Unified Form Language and shown how the language and its associated 
algorithms allow compact, readable and efficient specification of mathematical expressions related 
to variational formulations of partial differential equations. We have presented both high-level and 
detailed views of UFL to communicate its practical use and to provide developers and technical 
users a firm grounding in the design principles of UFL for understanding and building upon UFL. 

UFL is a stand-alone Python module that has been extensively used as part of the FEniCS soft- 
ware pipeline since 2009. The UFL functionality has been crucial in enabling advanced automated 
finite element algorithms in the FEniCS context, especially for complicated coupled systems of 
equations and for problems for which automatic differentiation dramatically reduced the burden on 
the application developer. UFL has also proven to be extensible beyond the core implementation, 
as exemplified by Nikbakht and Wells [351 an d Massing et al. [32] in the context of extended finite 
element methods, and Markall et al. [3U] in relation to code generation for different architectures. 
UFL is an actively developed project and continues to further extend the power and expressiveness 
of the language. 



Acknowledgments 

The authors wish to thank Kent- Andre Mardal and Johannes Ring for their contributions to UFL, 
and Pearu Peterson for discussions about symbolic representations during the initial design phase. 



References 

[1] Boost, 2012. URL |http : //www . boost . org/' 



[2] DOLFIN, 2012. URL http://launchpad.net/dolfin 



[3] PyChecker: A Python source code checking tool, 2012. URL http://pychecker. 



sourcef orge .net. 



[4] UFL, 2012. URL http://launchpad.net/ufl 



[5] M. S. Alnaes. A compiler framework for automatic linearization and efficient discretization 
of nonlinear partial differential equations. PhD thesis, University of Oslo, 2009. 

[6] M. S. Alnaes and K.-A. Mardal. On the efficiency of symbolic computations combined with 
code generation for finite element methods. ACM Trans. Math. Software, 37(1), 2010. 

[7] M. S. Alnaes and K.-A. Mardal. SyFi and SFC: Symbolic finite elements and form compilation. 
In A. Logg, K.-A. Mardal, and G. N. Wells, editors, Automated Solution of Differential Equa- 
tions by the Finite Element Method, volume 84 of Lecture Notes in Computational Science 
and Engineering, chapter 15. Springer, 2012. 

[8] M. S. Alnaes, A. Logg, K.-A. Mardal, O. Skavhaug, and H. P. Langtangen. Unified framework 
for finite element assembly. Int. J. Computational Science and Engineering, 4(4):231-244, 



2009. URL http : //dx . doi . org/10 . 1504/1 JCSE . 2009 . 029160 



[9] M. S. Alnass, A. Logg, and K.-A. Mardal. UFC: a finite element code generation interface. In 
A. Logg, K.-A. Mardal, and G. N. Wells, editors, Automated Solution of Differential Equations 
by the Finite Element Method, volume 84 of Lecture Notes in Computational Science and 
Engineering, chapter 16. Springer, 2012. 

[10] D. N. Arnold. An interior penalty finite element method with discontinuous elements. SIAM 
J. Numer. Anal, 19(4):742-760, 1982. 

[11] D. N. Arnold, F. Brczzi, and M. Fortin. A stable finite element for the Stokes equations. 
Calcolo, 21(4):337-344, 1984. URL |http : //dx . doi . org/ 10 . 1007/BF0257617l| 



32 



[12] D. N. Arnold, R. S. Falk, and R. Winther. Finite element exterior calculus, homological 
techniques, and applications. Acta Numer., 15:1-155, 2006. 

[13] B. Bagheri and L. R. Scott. About Analysa. Technical Report TR-2004-09, Departments 
of Computer Science, University of Chicago, 2004. URL http://www.cs.uchicago.edu/ 
research/publications/techreports/TR-2004-09 

[14] G. Baumgartner, A. Auer, D. E. Bernholdt, A. Bibireata, V. Choppella, D. Cociorva, X. Gao, 
R. J. Harrison, S. Hirata, S. Krishnamoorthy, S. Krishnan, C. Lam, Q. Lu, M. Nooijen, R. M. 
Pitzer, J. Ramanujam, P. Sadayappan, and A. Sibiryakov. Synthesis of high-performance 
parallel programs for a class of ab Initio quantum chemistry models. Proceedings of the 
IEEE, 93(2):276-292, 2005. URL |http : //dx . doi . org/10 . 11 09/ JPRO C . 2004 84031 1| 

[15] F. Brezzi, J. Douglas, and L. D. Marini. Two families of mixed elements for second order 
elliptic problems. Numer. Math., 47:217-235, 1985. 

[16] P. E. Farrcll, D. A. Ham, S. W. Funke, and M. E. Rognes. Automated derivation of the 
adjoint of high-level transient finite element programs. 2012. URL http://arxiv.org/abs/ 
11204.55771 

[17] E. Gamma, R. Helm, R. Johnson, and J. Vlissides. Design patterns: Abstraction and reuse 
of object-oriented design. ECOOP93-Object-Oriented Programming, pages 406-431, 1993. 

[18] A. Griewank. On automatic differentiation. In M. Iri and K. Tanabe, editors, Mathemat- 
ical Programming: Recent Developments and Applications, pages 83-108. Kluwer Academic 
Publishers, 1989. 

[19] J. Karczmarczuk. Functional differentiation of computer programs. Higher-Order 
and Symbolic Computation, 14(l):35-57, 2001. URL |http : //dx . doi . org/ 10 . 10237a7 
110115012321971 

[20] R. C. Kirby and A. Logg. A compiler for variational forms. A CM Trans. Math. Software, 32 
(3), 2006. URL |http : //dx . doi . org/10 . 1 145/1 163641 . 1163644| 

[21] R. C. Kirby and A. Logg. Finite element variational forms. In A. Logg, K.-A. Mardal, and 
G. N. Wells, editors, Automated Solution of Differential Equations by the Finite Element 
Method, volume 84 of Lecture Notes in Computational Science and Engineering, chapter 5. 
Springer, 2012. 

[22] J. Korelc. Automatic generation of finite-element code by simultaneous optimization of ex- 
pressions. Theor. Comput. Sci., 187(1 2):231 248, 1997. URL http : //dx . doi . org/10 . 1016/ 
IS0304-3975 (97) 00067-4| 

[23] R. J. Labeur and G. N. Wells. Energy stable and momentum conserving hybrid finite element 
method for the incompressible Navier-Stokes equations. SIAM J. Sci. Comput., 34(2):A889- 
A913, 2012. URL htt p : //dx . doi . org/10 . 1 137/100818583| 

[24] A. Logg and G. N. Wells. DOLFIN: Au tomated finite element computing. ACM Trans. Ma th. 
Software, 37(2) :20: 1-20:28, 2010. URL |http : //dx . doi . org/10 . 1145/1731022 . 1731030| 

[25] A. Logg, K.-A. Mardal, and G. N. Wells, editors. Automated Solution of Differential Equations 
by the Finite Element Method, volume 84 of Lecture Notes in Computational Science and 
Engineering. Springer, 2012. 

[26] A. Logg, K. B. Glgaard, M. E. Rognes, and G. N. Wells. FFC: The FEniCS form compiler. In 
A. Logg, K.-A. Mardal, and G. N. Wells, editors, Automated Solution of Differential Equations 
by the Finite Element Method, volume 84 of Lecture Notes in Computational Science and 
Engineering, chapter 11. Springer, 2012. 



33 



A. Logg, G. N. Wells, and J. Hake. DOLFIN: A C++/Python finite element library. In 
Automated Solution of Differential Equations by the Finite Element Method, volume 84 of 
Lecture Notes in Computational Science and Engineering, chapter 10. Springer, 2012. 

K. Long, R. Kirby, and B. Van Bloemen Waanders. Unified embedded parallel finite element 
computations via software-based Frechet differentiation. SIAM J. Sci. Comput., 32:3323- 
3351, 2010. 

M. Maraldi, G. N. Wells, and L. Molari. Phase field model for coupled displacive and diffusive 
microstructural processes under thermal loading. J. Mech. Phys. Solids, 59(8): 1596-1612, 
2011. URL |http : //dx . doi . org/10 . 1016/j . jmps . 2011 . 04 . Olfj 

G. Markall, A. Slemmer, D. Ham, P. Kelly, C. Cantwell, and S. Sherwin. Finite element 
assembly strategies on multi-core and many-core architectures. Internat. J. Numer. Methods 
Fluids, 2012. URL |http : //dx . doi . org/10 . 1002/f Id . 3648"! 

G. R. Markall, D. A. Ham, and P. H. J. Kelly. Towards generating optimised finite element 
solvers for GPUs from high-level specifications. Proc. Comp. Sci., 1(1):1815-1823, 2010. 

A. Massing, M. G. Larson, A. Logg, and M. E. Rognes. A stabilized Nitsche fictitious domain 
method for the Stokes problem. 2012. URL |http : //arxiv . org/abs/1206 . 1933 

M. Mortensen, H. P. Langtangen, and G. N. Wells. A FEniCS-bascd programming framework 
for modeling turbulent flow by the Reynolds-averaged Navier-Stokes equations. Adv. Water 
Res., 34(9):1082-1101, 2011. URL |http : //dx . doi . org/j advwatres . 2011 . 02 . 013| 

J. C. Nedelec. A new family of mixed finite elements in M. 3 . Numer. Math., 50(1):57-81, 1986. 

M. Nikbakht and G. N. Wells. Automated modelling of evolving discontinuities. Algorithms, 
2(3):1008-1030, 2009. URL |http : //dx . doi . org/10 . 3390/a2031008| 

K. B. 01gaard and G. N. Wells. Optimisations for quadrature representations of finite element 
tensors through automated code generation. ACM Trans. Math. Software, 37(l):8:l-8:23, 
2010. 

K. B. Olgaard, A. Logg, and G. N. Wells. Automated code generation for discontinuous 



Galerkin methods. SIAM J. Sci. Comput, 31(2):849-864, 2008. URL http : //dx . doi . org/ 
|10.1137/070710032i 

B. A. Pearlmutter and J. M. Siskind. Lazy multivariate higher-order forward-mode AD. 
In Proceedings of the 34th annual ACM SIGPLAN-SIGACT symposium on Principles of 
programming languages, pages 155-160, Nice, France, 2007. URL [http : / /dxTd oi . org / 10 . | 
11145/1190216. 11902421 

C. Prud'homme. A domain specific embedded language in C++ for automatic differentiation, 
projection, integration and variational formulations. Scientific Programming, 14:81-110, 2006. 



C. Prud'homme. Feel++, 2011. URL |https : //forge . imag . f r/pro j ects/lif e/| 



M. E. Rognes and A. Logg. Automated goal-oriented error control I: stationary variational 
problems. 2012. URL |http : //arxiv . org/abs/1204 . 6643[ 

M. E. Rognes, R. C. Kirby, and L. A. Efficient assembly of H(div) and H (curl) conforming 
finite elements. SIAM J. Sci. Comput, 31(6):4130-4151, 2009. URL http://dx .doi .org/| 
110.1137/ 08073901X1 

E. Rosseel and G. N. Wells. Optimal control with stochastic PDE constraints and uncertain 
controls. Comput Methods Appl. Mech. Engrg., 213-216:152-167, 2012. URL |http : //dx. | 



doi . org/10 . 1016/j . cma. 2011 . 11 . 026 



34 



[44] J. M. Siskind and B. A. Pearlmutter. Nesting forward-mode AD in a functional framework. 
Higher Order Symbol. Comput, 21(4):361-376, 2008. URL http : //dx . doi . org/10 . 1007/ 
S10990-008-9037-1. 



[45] C. Taylor and P. Hood. A numerical solution of the Navier-Stokes equations using the finite 
element technique. Comput. Fluids., 1(1):73-100, 1973. 

[46] P. Wang. FINGER: A symbolic system for automatic generation of numerical programs in 
finite element analysis. J. Symbolic Comput., 2(3):305-316, 1986. URL [http://dx.doi - org/ 
|10 . 1016/S0747-7171 (86)80029 :: 3l 

[47] G. N. Wells. Analysis of an interface stabilised finite element method: The advection-diffusion- 
reaction equation. SI AM J. Numer. Anal, 49(1):87-109, 2011. URL |http : //dx . doi . org/| 
110.1137/0907754641 

[48] J. Xiong, J. Johnson, R. Johnson, and D. Padua. SPL: A language and compiler for DSP 
algorithms. SIGPLAN Not., 36:298-308, 2001. URL |http : //dx . doi . org/ 10 . 1 145/381694 . 
378860. 



35 



